Complete identification of alkali sites in ion conducting lithium silicate glasses: a 

computer study of ion dynamics 
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The available sites for ions in a typical disordered ionic conductor are determined. For this 
purpose we devised a straightforward algorithm which via cluster analysis identifies these sites from 
long time ionic trajectories below the glass transition. This is exemplified for a lithium silicate glass 
(Li20) a .(Si02)( 1 _ :z .- ) for x = 0.5 and x = 0.1. We find for both concentrations that the number of 
sites is only slightly bigger than the number of ions. This result suggests a theoretical description 
of the dynamics in terms of mobile vacancies as most appropriate. Furthermore identification of 
the ionic sites allows one to obtain detailed characteristics of the ionic motion, e.g. quantification 
of correlated forward-backward jumps. 



Ion conducting glasses have been investigated by vari- 
ous experimental methods, including EXAFSEIQ, NMR 
[J and conductivity spectroscopy 0, 07 Whereas 
quite detailed information about the local structure has 
become available the mechanism of dynamics is still un- 
der debate although consensus has been reached that ion 
dynamics can be described as jumps of the mobile ions 
Hi 0- El • Whereas some authors stress the relevance of 
the disordered energy landscape 0, 0] supplied by the 
network, others relate the complexity of ion dynamics to 
the Coulomb interaction among mobile ions |fj . Further- 
more it has been argued from structural considerations 
that the distribution of alkali ions is inhomogeneous • 

For a closer understanding of ion dynamics microscopic 
information as supplied by molecular dynamics (MD) 
simulations is highly welcome. Jund et al. 13] found 
preferential pathways as a dynamical phenomenon by 
counting the number of different alkali ions that passed 
through subvolumes of the simulation box. The result- 
ing subvolumes visited by the largest number of different 
alkali ions may be interpreted as fast areas. They form a 
network of conduction paths. Oviedo and Sanz have ar- 
gued on a qualitative basis that for alkali concentrations 
lower than 10% the alkali ions are always surrounded by 
non-bridging oxy gens (NBOs) despite the overall small 
number of NBOs [14j . This implies that new NBOs are 
formed via breaking of Si-0 bonds along the alkali trajec- 
tory. In contrast, for higher concentrations hopping dy- 
namics between so-called micro-channels is proposed such 
that no formation of new NBOs would be necessary. Hor- 
bach |l5| showed in simulations of (Na20)2(Si02), that 
the coherent intermediate scattering function of sodium 
relaxes only on the same timescale as those of the net- 
work species, proving the existence of stable alkali sites 
below the glass transition. Cormack 0] et al. have re- 
cently investigated the mechanism of sodium migration 
in simulations of (Na2O) 2s(SiO2) 75 glasses, observing 
a few sequences of jumps between selected sites. He inter- 
preted the resulting dynamics as the motion of vacancies 



and pointed out that the identification of all sites in the 
glass, empty or populated at a given time, would be most 
useful for obtaining a deeper understanding of the kinetic 
mechanisms. 

Determination of all alkali sites in the silicate network 
is the purpose of this work. These sites can be defined as 
any regions in the network where alkali ions stay. They 
must correspond to local minima of the potential energy 
for the ions. Thus determination of the ionic sites is 
equivalent to determining the effective energy landscape 
experienced by the ions. The potential energy involves 
the network contribution as well as the averaged interac- 
tion with other ions. 

We analyse lithium silicate glasses (Li 2 0) :J ,(Si02)( 1 _ ;1 ,) 
with x — 0.5 and x = 0.1 which are common model sys- 
tems for computer investigations of ion dynamics. The 
system size has been chosen such that for x — 0.5 we 
have 384 lithium ions and for x — 0.1 80 lithium ions. 
Simulations are performed in the NVT ensemble at tem- 
peratures 640 K and 920 K and simulation box lengths 
22.96 A and 25.28 A, respectively. The resulting densi- 
ties match experimental values |l7j . The temperatures 
have been chosen such that the lithium diffusion con- 
stants in both systems have similar values. We use a 
pair potential developed by J. Habasaki [l^. Previous 
studies with this potential have shown good agreement 
with ex peri mental results for static and dynamic quanti- 
ties 19, 20, 21] . The molecular dynamics simulations are 
performed with a modified version of Moldy 22]. Since 
both temperatures are below the respective glass transi- 
tion temperatures we proceed in two steps. We start from 
a configuration which was first equilibrated just above 
the glass transition. It is then propagated at the selected 
temperature for ca. 20 ns in the NVT ensemble until the 
lithium subsystem is equilibrated. Afterwards we per- 
form a production run for approximately 50 ns. During 
this run, positions of all particles are recorded every 0.1 
ps. More details about the simulations can be found in 
\2l\ . We checked the validity of our results for indepen- 



0,1 



ul I I I 




Ij-S-B □ □ □ B-B-B-H-H-B-B B -B-H-EB 



10 100 1000 10000 
t(ps) 



FIG. 1: Mean square displacements of lithium and oxygen in 
the simulations. Included is a line with slope 1. 



dent data of up to 10 ns duration. This time is long 
enough for the lithium dynamics to become diffusive, so 
that just the dynamics relevant for ionic transport is in- 
cluded. 

In Fig. [I]we report the mean square displacement of 
oxygen and lithium. The dynamics of the lithium ions be- 
comes diffusive for tdiff ~ Ins for x = 0.5 and tdiff ~ 10 
ns for x = 0.1, as obtained from longer simulations. Evi- 
dently, the oxygen ions are localised on the ns-time scale. 
The network of silicon and oxygen atoms and with this 
the lithium sites can therefore be considered as static. 
As a consequence the ionic sites are basically unmodified 
during our run. Interestingly, as shown recently 0> US 
the remaining local fluctuations of the network are nev- 
ertheless essential in promoting the ionic jumps. 

Based on the lithium trajectories we have devised a 
straightforward algorithm to identify ionic sites. As a ba- 
sic idea we let the ions decide via their dynamics where 
good sites are rather than determining the potential en- 
ergy, as done e.g. in [23J. First the simulation box is 
divided into cubic cells with size approx. (0.3 A) 3 . They 
are small enough to enable resolution of the shape of 
the ionic sites. During our MD-run we count how many 
time steps a cell is visited by a lithium ion. The cells 
with nonzero counts describe the portion of the system 
that has been visited by lithium ions. These cells include 
the lithium sites as well as the connecting paths between 
them. To identify the ionic sites and to eliminate the 
paths between the sites, cells with less than a minimum 
value M of counts are dismissed. With the remaining 
cells a cluster analysis is performed. Cells which share 
a face are grouped into one cluster. Cells that are not 
directly or indirectly connected over common faces form 
different clusters. The value of M has been determined 
by the condition that the number of distinct ionic clusters 
be a maximum. For a smaller threshold some lower pop- 
ulated clusters would not have been detected, for a larger 




FIG. 2: Rendering of a slice of the simulation box for 
(Li a O) a .(Si0 2 ) 1 _ !1 . with x = 0.5(top) and x = O.l(bottom), 
together with clusters of cells with high average lithium oc- 
cupation. (Created with the software vmd |2(|.) 



threshold distinct clusters would be connected to single 
clusters because the transition path would be included. 
With M chosen according to this condition, 2.9% and 
1.3% of all cells are included for x = 0.5 and x = 0.1, 
which contain ss 90% of all lithium positions in both 
cases. The precise value of M is of no relevance for the 
latter analysis. 

For x = 0.5 and x = 0.1 respectively, 378 and 76 clus- 
ters are found. The result is illustrated in Fig. [3 showing 
a snapshot of a slice through each system. The silicate 
network is shown as tubes, with silicon bright and oxy- 
gens dark. Lithium ions are depicted as large spheres. 
The white objects are the clusters obtained from the 
analysis described above, with each small white sphere 
representing one of the cells. It can be seen that the clus- 
ters are in fact discrete objects of compact shape. Most of 
the lithium ions reside inside of one of these clusters, and 
most of the clusters are occupied by exactly one lithium 
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FIG. 3: Mean occupation of clusters. Inset: effective radius 
r e ff = (3V/47r) 1//3 of clusters, including the small clusters 
excluded from further analysis (see text). 



ion. As will be quantified further below some clusters 
can accommodate up to three lithium ions. 

For x=0.5 (x=0.1) in total 29 (14) additional small 
clusters have been excluded. They are set apart by a 
conspicious drop in the effective radius plot; see Fig. [3] 
Detailed analysis characterised most of them as simple 
satellites of bigger clusters. The remaining ones serve as 
short-lived saddle-like states during transitions between 
bigger clusters. The small clusters need thus not be con- 
sidered as sites themselves. 

On the basis of the clusters, the dynamics of a lithium 
ion can be expressed as successive residences in clusters 
separated by shorter time intervals where the lithium ion 
is outside of any cluster. When a lithium ion leaves one 
cluster and moves into a different one, this sequence of 
events is recorded as a jump. In contrast, if an ion ex- 
plores part of the volume outside a cluster without enter- 
ing a second one and comes back to the original cluster, 
it is regarded as occupying the same cluster during the 
whole time. Consequently, the residence time of a lithium 
ion in a cluster is defined as the time between its jump 
into this cluster and the subsequent jump to another clus- 
ter. With this definition on average a lithium ion in the 
system with x = 0.5 belongs to a cluster at 99.5% of all 
times. The remaining 0.5% correspond to jumps between 
different sites. For x — 0.1, the portion of times corre- 
sponding to jumps is 1.1%. The typical distance travelled 
by the lithium ions during a jump coincides well with the 
first maximum of the radial distribution function of the 
clusters in both cases. For x — 0.5, 93% of all jumps 
are inside the first neighbour shell (94% for x = 0.1) as 
defined by the first minimum of the radial distribution 
function. Thus the ion dynamics can be described as 
jumps between adjacent clusters. 

Having defined the residence times the mean number 
of lithium ions occupying each cluster can easily be eval- 



uated. It is plotted in Fig. [3] for x = 0.5, with the 
clusters ordered from left to right by decreasing volume 
V (the effective radius r e ff = (SV/Att) 1 ^ 3 of the individ- 
ual clusters is plotted in the inset). The occupation data 
confirms the impression gained from Fig. [21 For most 
clusters the value is close to unity. These clusters are 
occupied by a single ion during the largest part of the 
simulation. Their effective radius r e ff ~ 0.6A is very 
similar. Only the largest clusters offer space for more 
than one lithium ion, i.e. they contain more than one 
site, and often they are indeed multiply occupied. 

Generally an empty site must be available for a jump of 
a lithium ion. The total number of sites can now be esti- 
mated from the occupation data in Fig. |21 For the larger 
clusters with cluster index < 19 we have taken the next 
higher integer number to estimate the number of avail- 
able sites. In this way we find 3 clusters with 3 sites and 
9 clusters with 2 sites. Starting from cluster index 14 one 
finds a few clusters with occupation number 1 + e with 
e < 0.2. Here it would be rather unphysical to attribute 
two sites to this cluster. Rather a more detailed analy- 
sis of explicit ion trajectories has revealed the following 
scenario. At most time instances a single ion stayed in 
this cluster. During the short times during which this 
ion explored the immediate neighborhood of this cluster 
without entering a new cluster (typical excursion length 
scales are 2 A from the center of the cluster) a second 
ion may briefly enter this cluster and immediately leave 
this cluster again. Thus one may either say that these 
clusters only contain a single site or, to be more conser- 
vative with respect to finding an upper bound of the total 
number of sites, attribute 1 + e sites to this cluster. With 
the latter variant we obtain 395 sites for our 384 lithium 
ions. Analogously, 86 sites were found for the 80 ions for 
x = 0.1. These results shed new light on the mechanism 
of ion dynamics. Rather than speaking of individual ions 
jumping between accessible sites it is more appropriate 
to speak to first approximation of vacancy dynamics, i.e. 
the dynamics of the non-interacting unoccupied sites in 
the energy landscape. 

Actually, Dyre [23 has recently suggested this type 
of scenario based on general reasoning. He relates the 
possible near-equality of sites and ions to the effect that 
during network formation at the glass transition only a 
minimum number of sites will be formed due to ener- 
getic reasons. This scenario is also compatible with the 
counter-ion model, proposed by Dieterich and coworkers 
|2sj . but is at variance with single-particle approaches. 

In the remaining part of this paper we would like to 
characterize the ion dynamics in terms of the individ- 
ual clusters. Fig. ^ shows the cumulative distribution 
function S(t) of residence times. Thus S(t) denotes the 
fraction of clusters with an average residence time less 
than r. Its median T medlan , defined by S(T median ) = 0.5, 
is close to 1 ns for x = 0.5 and 250 ps for x = 0.1. These 
numbers have to be compared with the median of the 
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FIG. 4: Cumulative distribution function S(t) of residence 
times and time dependence of the backjump probability 
Pback{r) for systems with x = 0.5 and x — 0.1. Arrows mark 
the times at which the diffusive regime is reached. 



transition time for a jump between two clusters, which is 
500 fs for x = 0.5 and 200 fs for x = 0.1. This dramatic 
time scale separation clearly justifies the hopping picture 
in this kind of materials, as known from previous simu- 
lations |23. The lithium dynamics on long time scales 
can thus be mapped to jumps between the clusters. The 
generally higher frequency and lower duration of jumps 
at x = 0.1 can be traced to the higher temperature that 
was chosen for this system to achieve similar diffusivities. 

The effectiveness of jumps for diffusive transport is 
strongly decreased by correlated forward-backward mo- 
tion. We determined the probability pback that a jump 
from a cluster A to another cluster B is followed by a di- 
rect backjump to the cluster A. The result is also shown 
in Fig. 0] One can clearly see that clusters with short 
residence times show a stronger tendency for a corre- 
lated forward-backward jump. For the fastest clusters, 
Pback reaches more than 0.9, which is several times higher 
than the statistical value given by the inverse number 
of nearest neighbors (~ 0.2... 0.3). The difference in 
Pback(r m edian), which is w 0.8 for x = 0.1 compared to 
~ 0.55 for x = 0.5, may explain why the former system 
becomes diffusive only after most ions have performed 
several jumps, as shown by the ratio T me dian/tdif / >> 1. 
In contrast one finds T me cn an /tdiff ~ 1 for x — 0.5, mean- 
ing that only half of the ions have jumped at least once at 
the onset of diffusion. The dispersive regime at shorter 
times is thus dominated by the fastest lithium ions. 

Beyond the results shown in Fig. 0] a detailed anal- 
ysis of the dynamics, i. e. in terms of cooperativity, 
becomes feasible. Our approach may easily be gener- 
alised to more complex systems and important questions 
concerning, e.g., the character and the lifetime of the in- 
dividual sites may be analysed. 
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